#########################################################################################################
# Script Name: Power Analysis - Ethnic Mediation Paper
#
# Purpose: Run a simulation to perform a power analysis for the eliminated effect
#
# Author: Erica Ann Metheney
#
# Contact: data@gld.gu.se
#
# Notes: 
#          Y = stated vote decision
#          x1 = co-ethnic (randomly assigned)
#          x2 = natural mediator indicator  (randomly assigned)
#          r.x1.x2 = 0 (since both were independently and randomly assigned)
#          r.x1.y = -0.1  (by the theory should be negative. Exact value chosen from sample estimates)
#          r.x2.y = -0.1  (by the theory should be negative. Exact value chosen from sample estimates)
#
#########################################################################################################

# NECESSARY PACKAGES
library(InteractionPoweR)
library(ggplot2)
library(ggformula)

# GLOBAL VARIABLES
my.alpha = 0.05
my.r.x1.x2 = 0


# Combined Mediation
my.n = 7000

test_power_med<-power_interaction(
  n.iter = 10000,            # number of simulations
  alpha = my.alpha,          # alpha, for the power analysis
  N = my.n,                  # sample size
  r.x1x2.y = .05,            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,             # correlation between x1 and y
  r.x2.y = -0.1,             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,      # correlation between x1 and x2
  k.x1 = 2,                  # x1 is binary
  k.x2 = 2,                  # x2 is binary
  k.y = 2                    # y is binary
)


# Kenya Only Mediation
my.n = 800

test_power_Konly<-power_interaction(
  n.iter = 10000,            # number of simulations
  alpha = my.alpha,          # alpha, for the power analysis
  N = my.n,                  # sample size
  r.x1x2.y = .05,            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,             # correlation between x1 and y
  r.x2.y = -0.1,             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,      # correlation between x1 and x2
  k.x1 = 2,                  # x1 is binary
  k.x2 = 2,                  # x2 is binary
  k.y = 2                    # y is binary
)


# Malawi Only Mediation
my.n = 3900

test_power_Monly<-power_interaction(
  n.iter = 10000,            # number of simulations
  alpha = my.alpha,          # alpha, for the power analysis
  N = my.n,                  # sample size
  r.x1x2.y = .05,            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,             # correlation between x1 and y
  r.x2.y = -0.1,             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,      # correlation between x1 and x2
  k.x1 = 2,                  # x1 is binary
  k.x2 = 2,                  # x2 is binary
  k.y = 2                    # y is binary
)

# Zambia Only Mediation
my.n = 2200

test_power_Zonly<-power_interaction(
  n.iter = 10000,            # number of simulations
  alpha = my.alpha,          # alpha, for the power analysis
  N = my.n,                  # sample size
  r.x1x2.y = .05,            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,             # correlation between x1 and y
  r.x2.y = -0.1,             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,      # correlation between x1 and x2
  k.x1 = 2,                  # x1 is binary
  k.x2 = 2,                  # x2 is binary
  k.y = 2                    # y is binary
)


# Combined Results for Graph
test_power_full<-power_interaction(
  n.iter = 1000,                                            # number of simulations
  alpha = my.alpha,                                          # alpha, for the power analysis
  N = sort(unique(c(c(800,2200,3900,8000),floor(seq(from = 500, to = 10000, length.out = 20))))),    # sample size
  r.x1x2.y = .05,                                            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,                                             # correlation between x1 and y
  r.x2.y = -0.1,                                             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,                                      # correlation between x1 and x2
  k.x1 = 2,                                                  # x1 is binary
  k.x2 = 2,                                                  # x2 is binary
  k.y = 2                                                    # y is binary
)

#----Plot
test_power_full$pwr = 100*test_power_full$pwr
n <- 20
df = data.frame(x1 = c(800,2200,3900,8000),y1 = c(0,0,0,0),x2 = c(800,2200,3900,8000) , y2 = test_power_full$pwr[test_power_full$N %in% c(800,2200,3900,8000)])
df2 = data.frame(x1 = c(0,0,0,0),y1 =  test_power_full$pwr[test_power_full$N %in% c(800,2200,3900,8000)],x2 = c(800,2200,3900,8000) , y2 =  test_power_full$pwr[test_power_full$N %in% c(800,2200,3900,8000)])
ggplot_formatted <- ggplot(test_power_full) + 
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 1,color = "red", data = df)+ #vertical lines from sample size to powercurve
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 1,color = "red", data = df2)+ #horizontal lines from powercurve to y-axis
  geom_line(aes(x = N,y = pwr),size = 2)+
  ggtitle("Power Analysis for Detecting an Eliminated Effect of 0.05\n")+
  xlab("\nSample Size")+
  ylab("Power\n")+
  scale_x_continuous(limits = c(0,9000),expand = c(0,0),breaks=c(800,2200,3900,8000),labels=c("Kenya\n(n = 800)","Zambia\n(n = 2200)","Malawi\n(n = 3900)","Combined\n(n = 8000)"))+
  scale_y_continuous(limits = c(0,100),expand = c(0,0),breaks = test_power_full$pwr[test_power_full$N %in% c(800,2200,3900,8000)],labels=paste(test_power_full$pwr[test_power_full$N %in% c(800,2200,3900,8000)],"%",sep = ""))+
  theme_classic()+
  theme(axis.text.y = element_text(face = "bold",size = 8),
        axis.text.x = element_text(face = "bold",size = 8),
        axis.title.x = element_text(face = "bold",size = 10),
        axis.title.y = element_text(face = "bold",size = 10),
        plot.title = element_text(face = "bold",size = 12))

ggsave(filename = "PowerAnalysis_EliminatedEffect.png",
       plot = ggplot_formatted,
       device = "png",
       width = 6,
       units = "in",
       dpi = 600)





test_power_full2 <-power_interaction(
  n.iter = 1000,                                            # number of simulations
  alpha = my.alpha,                                          # alpha, for the power analysis
  N = sort(unique(c(c(800,2200,3900,8000),floor(seq(from = 500, to = 10000, length.out = 20))))),    # sample size
  r.x1x2.y = c(0.02,.05),                                            # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = -0.1,                                             # correlation between x1 and y
  r.x2.y = -0.1,                                             # correlation between x2 and y
  r.x1.x2 = my.r.x1.x2,                                      # correlation between x1 and x2
  k.x1 = 2,                                                  # x1 is binary
  k.x2 = 2,                                                  # x2 is binary
  k.y = 2                                                    # y is binary
)


test_power_full2$pwr = 100*test_power_full2$pwr
test_power_full2$`Effect Size` = factor(as.character(test_power_full2$r.x1x2.y), levels = c("0.05","0.02"))
test_power_full2.02 = subset(test_power_full2,`Effect Size` == "0.02")
test_power_full2.05 = subset(test_power_full2,`Effect Size` == "0.05")
n <- 20
df.02 = data.frame(x1 = c(800,2200,3900,8000),y1 = c(0,0,0,0),x2 = c(800,2200,3900,8000) , y2 = test_power_full2.02$pwr[test_power_full2.02$N %in% c(800,2200,3900,8000)])
df2.02 = data.frame(x1 = c(0,0,0,0),y1 =  test_power_full2.02$pwr[test_power_full2.02$N %in% c(800,2200,3900,8000)],x2 = c(800,2200,3900,8000) , y2 =  test_power_full2.02$pwr[test_power_full2.02$N %in% c(800,2200,3900,8000)])
df.05 = data.frame(x1 = c(800,2200,3900,8000),y1 = c(0,0,0,0),x2 = c(800,2200,3900,8000) , y2 = test_power_full2.05$pwr[test_power_full2.05$N %in% c(800,2200,3900,8000)])
df2.05 = data.frame(x1 = c(0,0,0,0),y1 =  test_power_full2.05$pwr[test_power_full2.05$N %in% c(800,2200,3900,8000)],x2 = c(800,2200,3900,8000) , y2 =  test_power_full2.05$pwr[test_power_full2.05$N %in% c(800,2200,3900,8000)])

ggplot_formatted <- ggplot(test_power_full2) + 
  geom_line(aes(x = N,y = pwr,group = `Effect Size`,color = `Effect Size`,linetype = `Effect Size`),size = 1)+
  scale_color_manual(values = c("#EF9307","#79C9C6"))+
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 0.5,color = "#161B39", data = df.02)+ #vertical lines from sample size to powercurve
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 0.5,color = "#161B39", data = df2.02)+ #horizontal lines from powercurve to y-axis
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 0.5,color = "#161B39", data = df.05)+ #vertical lines from sample size to powercurve
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),linetype = 2,size = 0.5,color = "#161B39", data = df2.05)+ #horizontal lines from powercurve to y-axis
  ggtitle("Power Analysis for Detecting an Eliminated Effect\n")+
  xlab("\nSample Size")+
  ylab("Power\n")+
  scale_x_continuous(limits = c(0,9000),expand = c(0,0),breaks=c(800,2200,3900,8000),labels=c("Kenya\n(n = 800)","Zambia\n(n = 2200)","Malawi\n(n = 3900)","Combined\n(n = 8000)"))+
  scale_y_continuous(limits = c(0,100),expand = c(0,0),breaks = test_power_full2$pwr[test_power_full2$N %in% c(800,2200,3900,8000)],labels=paste(test_power_full2$pwr[test_power_full2$N %in% c(800,2200,3900,8000)],"%",sep = ""))+
  theme_classic()+
  theme(axis.text.y = element_text(face = "bold",size = 8),
        axis.text.x = element_text(face = "bold",size = 8),
        axis.title.x = element_text(face = "bold",size = 10),
        axis.title.y = element_text(face = "bold",size = 10),
        plot.title = element_text(face = "bold",size = 12))

ggsave(filename = "PowerAnalysis_EliminatedEffect.png",
       plot = ggplot_formatted,
       device = "png",
       width = 6,
       units = "in",
       dpi = 600)

